Question 1

#path.root <- "D:/OneDrive - University of Vermont/Classes/Spring2022/sdhmR/sdhmR-V2022.1" #Reed Laptop Path
path.root <- "C:/Users/14842/Documents/SDHM/sdhmR-V2022.1" #Lindsey Path
path.mod <- paste(path.root, "/data/exercise/traindat", sep = "")
path.preds <- paste(path.root, '/data/exercise/preds', sep = '')
path.maps <- paste(path.root, '/data/exercise/maps', sep = '')
path.gis <- paste(path.root, "/data/gis_layers", sep = "")
path.figs <- paste(path.root, "/powerpoints/figures", sep = "")

load libraries

  library(raster)
 options("rgdal_show_exportToProj4_warnings" = "none") # run this before library(rgdal)
  library(rgdal)
  library(gam)
  library(randomForest)
  library(dismo)
  library(gbm)
  library(sf)

START INITIALIZATION OF DATA STRUCTURES

load SDM models. I renamed the models and the threshold cuts to make it more clear raster stack of predictors assumes rasters all have same projection, extent & resolution

  setwd(path.preds)
  pers.list <- list.files(pattern = ".img$") # list of .img files; $ strips extra
  pers.list # examine
##  [1] "etpt_5.img"     "etpt_6.img"     "etpt_sprin.img" "exp1nrm.img"   
##  [5] "exp3nrm.img"    "exp5nrm.img"    "mind_yr_av.img" "prad_sw_di.img"
##  [9] "prec_w_hal.img" "prec_winte.img" "rough_1k.img"   "tave_s_hal.img"
## [13] "tave_sprin.img" "tmax_s_hal.img" "tmax_summe.img" "topos.img"
  pers.dom <- stack(pers.list) # build raster stack
  pers.dom # examine stack
## class      : RasterStack 
## dimensions : 1518, 1478, 2243604, 16  (nrow, ncol, ncell, nlayers)
## resolution : 0.008333333, 0.008333333  (x, y)
## extent     : -114.6167, -102.3, 29.49167, 42.14167  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## names      :       etpt_5,       etpt_6,   etpt_sprin,      exp1nrm,      exp3nrm,      exp5nrm,   mind_yr_av,   prad_sw_di,   prec_w_hal,   prec_winte,     rough_1k,   tave_s_hal,   tave_sprin,   tmax_s_hal,   tmax_summe, ... 
## min values :     0.000000,    69.000000,     0.000000, -1152.000000,  -812.000000,  -715.000000, -2294.333252,  6411.166016,     3.666667,     3.333333,     0.000000,    16.166666,   -86.000000,    63.833332,   100.000000, ... 
## max values :    100.00000,    104.00000,     90.33334,   1246.00000,    945.00000,    895.00000,    249.33333,  14441.50000,    136.83333,    140.33333,    243.85391,    308.33334,    226.66667,    392.33334,    424.00000, ...
  names(pers.dom) 
##  [1] "etpt_5"     "etpt_6"     "etpt_sprin" "exp1nrm"    "exp3nrm"   
##  [6] "exp5nrm"    "mind_yr_av" "prad_sw_di" "prec_w_hal" "prec_winte"
## [11] "rough_1k"   "tave_s_hal" "tave_sprin" "tmax_s_hal" "tmax_summe"
## [16] "topos"
# load SDM models. I renamed the models and the threshold cuts to make it more clear
 
setwd(paste(path.mod, sep = ""))
load('ex7.RData')
  LR.cut <- mod.cut$glm.cvpred ### Assuming this is the cut point? In his example it's just one value
  LR <- mod2.LR
load('ex8.RData')
  GAM.cut <- cut$pred
  GAM <- mod1.GAM
load('ex9.RData')
  MAX <- mod1.MAX
  MAX.cut <- mod.cut$spec_sens
load('ex10.RData')
  RF <- pers.RF
  RF.cut <- mod.cut$persRF.pred
load('ex11.RData')
  BRT <- pers.BRT
  BRT.cut <- mod.cut$pred
  

### Building a list of cut points.
  cut.list <- c(LR.cut,
                GAM.cut,
                MAX.cut,
                RF.cut,
                BRT.cut
  )
  cut.list # threshold cuts?
## [1] 0.5077240 0.6200000 0.4894119 0.3400000 0.5077240
# raster stack of predictors assumes rasters all have same projection, extent & resolution
  setwd(path.preds)
  pers.list <- list.files(pattern = ".img$") # list of .img files; $ strips extra
  pers.list # examine
##  [1] "etpt_5.img"     "etpt_6.img"     "etpt_sprin.img" "exp1nrm.img"   
##  [5] "exp3nrm.img"    "exp5nrm.img"    "mind_yr_av.img" "prad_sw_di.img"
##  [9] "prec_w_hal.img" "prec_winte.img" "rough_1k.img"   "tave_s_hal.img"
## [13] "tave_sprin.img" "tmax_s_hal.img" "tmax_summe.img" "topos.img"
  pers.dom <- stack(pers.list) # build raster stack
  pers.dom # examine stack
## class      : RasterStack 
## dimensions : 1518, 1478, 2243604, 16  (nrow, ncol, ncell, nlayers)
## resolution : 0.008333333, 0.008333333  (x, y)
## extent     : -114.6167, -102.3, 29.49167, 42.14167  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## names      :       etpt_5,       etpt_6,   etpt_sprin,      exp1nrm,      exp3nrm,      exp5nrm,   mind_yr_av,   prad_sw_di,   prec_w_hal,   prec_winte,     rough_1k,   tave_s_hal,   tave_sprin,   tmax_s_hal,   tmax_summe, ... 
## min values :     0.000000,    69.000000,     0.000000, -1152.000000,  -812.000000,  -715.000000, -2294.333252,  6411.166016,     3.666667,     3.333333,     0.000000,    16.166666,   -86.000000,    63.833332,   100.000000, ... 
## max values :    100.00000,    104.00000,     90.33334,   1246.00000,    945.00000,    895.00000,    249.33333,  14441.50000,    136.83333,    140.33333,    243.85391,    308.33334,    226.66667,    392.33334,    424.00000, ...
  names(pers.dom)
##  [1] "etpt_5"     "etpt_6"     "etpt_sprin" "exp1nrm"    "exp3nrm"   
##  [6] "exp5nrm"    "mind_yr_av" "prad_sw_di" "prec_w_hal" "prec_winte"
## [11] "rough_1k"   "tave_s_hal" "tave_sprin" "tmax_s_hal" "tmax_summe"
## [16] "topos"
# LR prediction and classification
  setwd(paste(path.maps,
              sep = ""))
  names(pers.dom) <- paste(names(pers.dom),'img')
  modFprob.LR <- predict(pers.dom,
                         LR,
                         filename = "modFprob.GLM.img", 
                         type = "response",
                         fun = predict,
                         index = 2,
                         overwrite = T)
  modFclass.LR <- reclassify(pers.prob,
                             filename = "modFclass.LR.img",
                           c(0,
                             mod.cut[[2]],
                             0,
                                        mod.cut[[2]],
                             1,
                             1),
                           overwrite=TRUE)
# giggle plots
  par(mfrow = c(1, 2))
  plot(modFprob.LR,
       axes = T,
       main = "LR probability map") # plot probability map
 
  plot(modFclass.LR,
       axes = T,
       main = "LR classfied map") # plot classified map

  par(mfrow = c(1, 1))
# GAM prediction and classification
  library(gam)
  modFprob.GAM <- predict(pers.dom,
                          GAM,
                          filename = "modFprob.GAM.img", 
                          type = "response",
                          fun = predict,
                          index = 2,
                          overwrite = T)
  modFclas.GAM <- reclassify(modFprob.GAM,
                             filename = "modFclas.GAM.img", 
                             (c(0,
                                GAM.cut,
                                0,
                                GAM.cut,
                                1,
                                1)),
                             overwrite = T) 
  # giggle plots
  par(mfrow = c(1, 2))
  plot(modFprob.GAM,
       axes = T,
       main = "GAM probability map") # plot probability map
  plot(modFclas.GAM,
       axes = T,
       main = "GAM classfied map") # plot classified map

  par(mfrow = c(1, 1)) 
# RF prediction and classification
  library(randomForest)  # load library
  modFprob.RF <- predict(pers.dom,
                         RF,
                         filename = "modFprob.RF.img", 
                         type = "prob",
                         fun = predict,
                         index = 2,
                         overwrite = T) # prob map
  modFclas.RF <- reclassify(modFprob.RF,
                            filename = "modFclas.RF.img", 
                            (c(0,
                               RF.cut,
                               0,
                               RF.cut,
                               1,
                               1)),
                            overwrite = T) # class map
  # giggle maps
  par(mfrow = c(1, 2))
  plot(modFprob.RF,
       axes = T,
       main = "RF probability map") # plot probability map
  plot(modFclas.RF,
       axes = T,
       main = "RF classfied map") # plot classified map

  par(mfrow = c(1, 1))  
library(dismo)
library(gbm)
modFprob.BRT <- predict(pers.dom,
                        BRT,
                        n.trees = BRT$gbm.call$best.trees, 
                          type = "response",
                        filename = "modFprob.BRT.img",
                        overwrite = T) # prob map
  modFclas.BRT <- reclassify(modFprob.BRT,
                             c(0, BRT.cut,
                               0, BRT.cut,
                               1,
                               1), 
                             filename = "modFclas.BRT.img",
                             overwrite = T) # clas map
  # giggle plots
  par(mfrow = c(1, 2))
  plot(modFprob.BRT,
       axes = T,
       main = "BRT probability map") # plot probability map
  plot(modFclas.BRT, axes = T,
       main = "BRT classfied map") # plot classified map

  par(mfrow = c(1, 2))
  
  # save BRT plots if desired 
  setwd(path.figs)
  savePlot(filename = "mod06fig04.pdf",
           type = "pdf")
## Error in savePlot(filename = "mod06fig04.pdf", type = "pdf"): can only copy from 'windows' devices
  ## MAXENT prediction and classification
  setwd(path.maps)
  library(dismo) # load library
  modFprob.MAX <- predict(MAX,
                          pers.dom,
                          filename = "modFprob.MAX.img",
                          overwrite = T)  
  modFclas.MAX <- reclassify(modFprob.MAX,
                             filename = "modFclas.MAX.img", 
                             c(0,
                               MAX.cut,
                               0,
                               MAX.cut,
                               1,
                               1),
                             overwrite = T) # clas map
  
  # giggle plots
  par(mfrow = c(1, 2))
  plot(modFprob.MAX, axes = T, main = "MAX probability map") # plot probability map

  plot(modFclas.MAX, axes = T, main = "MAX classfied map") # plot classified map

  par(mfrow = c(1, 1))
  
  # save MAX plots if desired 
  setwd(path.figs)
  savePlot(filename = "mod06fig05.pdf", type = "pdf")
## Error in savePlot(filename = "mod06fig05.pdf", type = "pdf"): can only copy from 'windows' devices

END PREDICTION PROBABILITY AND CLASSIFED MAPS START ENSEMBLE PROCESS # load probability & classified maps; create stacks

  setwd(paste(path.maps, sep = ""))
  pr.list <- unlist(unique(strsplit(list.files(pattern = "Fprob"),
                                    ".aux.xml")))
  pr.list # examine
## [1] "modFprob.BRT.img" "modFprob.GAM.img" "modFprob.GLM.img" "modFprob.MAX.img"
## [5] "modFprob.RF.img"
  prob.dom <- stack(pr.list) # raster stack prob maps
  prob.dom # examine
## class      : RasterStack 
## dimensions : 1518, 1478, 2243604, 5  (nrow, ncol, ncell, nlayers)
## resolution : 0.008333333, 0.008333333  (x, y)
## extent     : -114.6167, -102.3, 29.49167, 42.14167  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## names      : modFprob.BRT, modFprob.GAM, modFprob.GLM, modFprob.MAX,  modFprob.RF 
## min values : 2.186834e-01, 2.220446e-16, 2.220446e-16, 1.048621e-05, 0.000000e+00 
## max values :    0.6525138,    0.9871084,    0.9911523,    0.7762663,    1.0000000
  # cl.list <- list.files(pattern = "Fclas")
  cl.list <- unlist(unique(strsplit(list.files(pattern = "Fclas"),
                                    ".aux.xml")))
  cl.list # examine
## [1] "modFclas.BRT.img" "modFclas.GAM.img" "modFclas.LR.img"  "modFclas.MAX.img"
## [5] "modFclas.RF.img"  "modFclass.LR.img"
  clas.dom <- stack(cl.list) # raster stack classified maps
  clas.dom # examine
## class      : RasterStack 
## dimensions : 1518, 1478, 2243604, 6  (nrow, ncol, ncell, nlayers)
## resolution : 0.008333333, 0.008333333  (x, y)
## extent     : -114.6167, -102.3, 29.49167, 42.14167  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## names      : modFclas.BRT, modFclas.GAM, modFclas.LR, modFclas.MAX, modFclas.RF, modFclass.LR 
## min values :            0,            0,           0,            0,           0,            0 
## max values :            1,            1,           1,            1,           1,            1
  # standardize all prediction maps 0-1.0
  maxValue(prob.dom) # extract max value for each prob map
## [1] 0.6525138 0.9871084 0.9911523 0.7762663 1.0000000
  layers <- {} # initialize (empty) list of raster layers
  for (i in 1:length(names(prob.dom))) {
    m1 <- prob.dom[[i]] # get a prob map
    m2 <- 1/maxValue(prob.dom[[i]]) * prob.dom[[i]] # standardize all probs to max=1
    m3 <- unlist(strsplit(names(prob.dom[[i]]), "[.]")) # split prob layer name apart
    names(m2) <- paste(m3[1], "STD.", m3[2], sep = "")  # assign name to raster value
    assign(paste(m3[1], "STD.", m3[2], sep = ""), m2) # assign new name to standardized layer
    layers <- c(layers, get(paste(m3[1], "STD.", m3[2], sep = "")))
  }
  probSTD.dom <- stack(layers)
  maxValue(probSTD.dom) # extract max value for each prob map; standardized ?
## [1] 1 1 1 1 1
  # save stacks
  setwd(path.mod)
  save(prob.dom,
       probSTD.dom,
       clas.dom,
       file = "ensemble.dom.RData")
  ######################################################
  # descriptive stats on raster maps: mean & sum
  prob.mean <- mean(prob.dom) # mean prob map
  # access min max values of descriptive stat raster
  maxValue(prob.mean) # max pr.mean
## [1] 0.8468679
  minValue(prob.mean) # min pr.mean
## [1] 0.04380246
  probSTD.mean <- mean(probSTD.dom) # mean standardized prob map
  clas.sum <- sum(clas.dom) # sum of models by cell
  
  # giggle plots: mean & sum
  par(mfrow = c(1, 2))
  plot(prob.mean,
       axes = T,
       main = "MEAN probability map") # plot probability map
  plot(st_geometry(pegr6.mahog),
       add = T)
## Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'plot': object 'pegr6.mahog' not found
  plot(st_geometry(pegr6.frame),
       add = T) # make pretty
## Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'plot': object 'pegr6.frame' not found
  plot(clas.sum,
       axes = T,
       main = "Concordance CLASS map: Ramp") # plot concordance map

  plot(st_geometry(pegr6.mahog),
       add = T)
## Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'plot': object 'pegr6.mahog' not found
  plot(st_geometry(pegr6.frame),
       add = T) # make pretty 
## Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'plot': object 'pegr6.frame' not found
  par(mfrow = c(1, 1))
  
  # save plots if desired 
  setwd(path.figs)
  savePlot(filename = "mod06fig06.pdf", type = "pdf")
## Error in savePlot(filename = "mod06fig06.pdf", type = "pdf"): can only copy from 'windows' devices
  ####
  # arithmetic on raster maps: multiplication & subsetting
  prob.mean1 <- (clas.sum > 0) * prob.mean  # mean where class=1
  clas.sum1 <- clas.sum > 0 # sum where No. models >=1
  
  # giggle plots: means X classifications
  par(mfrow = c(1, 2))
  plot(prob.mean1,
       axes = T,
       main = "MEAN probability map: Presence=1") # plot probability map
  plot(st_geometry(pegr6.mahog),
       add = T)
## Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'plot': object 'pegr6.mahog' not found
  plot(st_geometry(pegr6.frame),
       add = T) # make pretty
## Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'plot': object 'pegr6.frame' not found
  plot(clas.sum1, axes = T,
       main = "Concordance CLASS map: Union") # plot concordance map

  plot(st_geometry(pegr6.mahog), add = T)
## Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'plot': object 'pegr6.mahog' not found
  plot(st_geometry(pegr6.frame), add = T) # make pretty 
## Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'plot': object 'pegr6.frame' not found
  par(mfrow = c(1, 1))
  
  
  ####
  # calculate measures of variation
  prob.sd <- calc(prob.dom, sd) # sd prob map
  prob.var <- calc(prob.dom, var) # var prob map
  prob.cv <- (prob.sd/prob.mean)*100 # coefficient of variation
  
  # giggle plots: sd and var
  par(mfrow = c(1, 2))
  plot(prob.sd,
       axes = T,
       main = "SD probability map") # plot classified map
  plot(st_geometry(pegr6.mahog), add = T)
## Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'plot': object 'pegr6.mahog' not found
  plot(st_geometry(pegr6.frame), add = T) # make pretty
## Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'plot': object 'pegr6.frame' not found
  plot(prob.var, axes = T, main = "VAR probability map") # plot classified map

  plot(st_geometry(pegr6.mahog), add = T)
## Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'plot': object 'pegr6.mahog' not found
  plot(st_geometry(pegr6.frame), add = T) # make pretty 
## Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'plot': object 'pegr6.frame' not found

# make pretty

# save plots if desired #setwd(path.figs) #savePlot(filename = “mod06fig13.pdf”, type = “pdf”) ######## END ENSEMBLE INTERPRETATION ################################################################################ ########################################################################### ## Question #2 # Calculate the frequencies of “presence” in each of the 5 SDHMs

#### # some summary statistics: frequencies

  clas.freqM <- freq(clas.dom) # [0,1] freqs by clas map; rtns list
  clas.freqM[c(1,5)] # examine freqs for 1st 2 models
## $modFclas.BRT
##      value   count
## [1,]     0 1326676
## [2,]     1  865593
## [3,]    NA   51335
## 
## $modFclas.RF
##      value   count
## [1,]     0 1409045
## [2,]     1  783224
## [3,]    NA   51335
  clas.freqS <- data.frame(freq(clas.sum)) # freqs of concordance of models
  clas.freqS # examine; value is concordance freq
##   value   count
## 1     0 1240510
## 2     1   74298
## 3     2  111049
## 4     3  115015
## 5     4   39812
## 6     5  144216
## 7     6  467369
## 8    NA   51335
  clas.freqS$count[6]/sum(clas.freqS$count[2:6]) # prop. models concordance=5 
## [1] 0.297727

######## START ENSEMBLE INTERPRETATION #### # 6 prob map plots

  par(mfrow = c(2, 3))
  plot(modFprob.LR,
       axes= T,
       main = "LR probability map") # LR prob map
  plot(st_geometry(pegr6.mahog),
       add = T)
## Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'plot': object 'pegr6.mahog' not found
  plot(st_geometry(pegr6.frame),
       add = T) # make pretty
## Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'plot': object 'pegr6.frame' not found
  plot(modFprob.GAM,
       axes = T,
       main = "GAM probability map") # GAM prob map
  plot(st_geometry(pegr6.mahog), add = T)
## Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'plot': object 'pegr6.mahog' not found
  plot(st_geometry(pegr6.frame), add = T) # make pretty
## Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'plot': object 'pegr6.frame' not found
  plot(modFprob.MAX,
       axes = T,
       main = "MAX probability map") # MAX prob map
  plot(st_geometry(pegr6.mahog), add = T)
## Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'plot': object 'pegr6.mahog' not found
  plot(st_geometry(pegr6.frame), add = T) # make pretty
## Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'plot': object 'pegr6.frame' not found
  plot(modFprob.RF,
       axes = T,
       main = "RF probability map") # RF prob map
  plot(st_geometry(pegr6.mahog), add = T)
## Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'plot': object 'pegr6.mahog' not found
  plot(st_geometry(pegr6.frame), add = T) # make pretty
## Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'plot': object 'pegr6.frame' not found
  plot(modFprob.BRT,
       axes = T,
       main = "BRT probability map") # BRT prob map
  plot(st_geometry(pegr6.mahog), add = T)
## Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'plot': object 'pegr6.mahog' not found
  plot(st_geometry(pegr6.frame), add = T) # make pretty 
## Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'plot': object 'pegr6.frame' not found
  plot(prob.mean1,
       axes = T,
       main = "MEAN probability map: \nPresence=1") # MEAN prob map

  plot(st_geometry(pegr6.mahog), add = T)
## Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'plot': object 'pegr6.mahog' not found
  plot(st_geometry(pegr6.frame), add = T) # make pretty 
## Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'plot': object 'pegr6.frame' not found

#### # 6 classified maps

  par(mfrow = c(2, 3))
  plot(modFclass.LR,
       axes = T,
       main = "LR classfied map") # LR classified map
  
  plot(modFclas.GAM, axes = T, main = "GAM classfied map") # GAM classified map
  
  plot(modFclas.MAX, axes = T, main = "MAX classfied map") # MAX classified map
  
  plot(modFclas.RF, axes = T, main = "RF classfied map") # RF classified map
  
  plot(modFclas.BRT, axes = T, main = "BRT classfied map") # BRT classified map
  
  plot(clas.sum, axes = T, main = "Concordance CLASS \nmap: Ramp") # plot concordance map

  # save plots if desired 
  #setwd(path.figs)
  #savePlot(filename = "mod06fig11.pdf", type = "pdf")
  ####
  # 6 classification union maps
  clas.sumU <- clas.sum > 0 # sum where No. models =union
  clas.sum1 <- clas.sum > 1 # sum where No. models =2+
  clas.sum2 <- clas.sum > 2 # sum where No. models =3+
  clas.sum3 <- clas.sum > 3 # sum where No. models =4+
  clas.sum4 <- clas.sum > 4 # sum where No. models =5 
  
  # giggle plots
  par(mfrow = c(2, 3))
  plot(clas.sum, axes = T, main = "Concordance CLASS \nmap: Ramp")  # RAMP concordance map
  
  plot(clas.sumU, axes = T, main = "Concordance CLASS \nmap: Union")  # UNION concordance map
  
  plot(clas.sum1, axes = T, main = "Concordance CLASS \nmap: 2+")  # 2+ concordance map
  
  plot(clas.sum2, axes = T, main = "Concordance CLASS \nmap: 3+")  # 3+ concordance map
  
  plot(clas.sum3, axes = T, main = "Concordance CLASS \nmap: 4+")  # 4+ concordance map
  
  plot(clas.sum4, axes = T, main = "Concordance CLASS \nmap: 5") # 5 concordance map

# save plots if desired #setwd(path.figs) #savePlot(filename = “mod06fig12.pdf”, type = “pdf”)

  ####
  # some variance analyses
  q1 <- summary(prob.sd)[2] # 1st quartile
  q2 <- summary(prob.sd)[3] # 2nd quartile
  q3 <- summary(prob.sd)[4] # 3rd quartile
  prob.sd1q <- (prob.sd >= q1) * clas.sum1 # where 2+ models agree
  prob.sd2q <- (prob.sd >= q2) * clas.sum1 # where 2+ models agree
  prob.sd3q <- (prob.sd >= q3) * clas.sum1 # where 2+ models agree
  
  # giggle plots
  par(mfrow = c(2, 2))
  plot(prob.sd, axes = T, main = "SD probability map") # SD prob map
  plot(st_geometry(pegr6.mahog), add = T)
  plot(st_geometry(pegr6.frame), add = T) # make pretty
  
  plot(prob.sd1q, axes = Tmain = "1st quartile SD map") # 1stQ prob map
  plot(st_geometry(pegr6.mahog), add = T)
  plot(st_geometry(pegr6.frame), add = T) # make pretty
  
  plot(prob.sd2q, axes = T, main = "2nd quartile SD map") # 2ndQ prob map
  plot(st_geometry(pegr6.mahog), add = T)
  plot(st_geometry(pegr6.frame), add = T) # make pretty
  
  plot(prob.sd3q, axes = T, main = "3rd quartile SD map") # 3rdQ prob map
  plot(st_geometry(pegr6.mahog), add = T)
  plot(st_geometry(pegr6.frame), add = T)
## Error: <text>:16:32: unexpected '='
## 15:   
## 16:   plot(prob.sd1q, axes = Tmain =
##                                    ^

Question #3

#Tally the frequencies of concordance after “clipping” and compare with above

#### clipping...
  # Need to clip pers.DOM 

clas.freqS <- data.frame(freq(clas.sum)) # freqs of concordance of models
clas.freqS # examine; value is concordance freq
##   value   count
## 1     0 1240510
## 2     1   74298
## 3     2  111049
## 4     3  115015
## 5     4   39812
## 6     5  144216
## 7     6  467369
## 8    NA   51335
clas.freqS$count[6]/sum(clas.freqS$count[2:6])
## [1] 0.297727
 # 6 classification union maps
  clas.sumU.freq <- data.frame(freq(clas.sumU)) # sum where No. models =union
  clas.sum1.freq <- data.frame(freq(clas.sum1)) # sum where No. models =2+
  clas.sum2.freq <- data.frame(freq(clas.sum2)) # sum where No. models =3+
  clas.sum3.freq <- data.frame(freq(clas.sum3)) # sum where No. models =4+
  clas.sum4.freq <- data.frame(freq(clas.sum4))# sum where No. models =5 

  
  
clas.sumU.freq$count
## [1] 1240510  951759   51335
  clas.sum1.freq$count 
## [1] 1314808  877461   51335
  clas.sum2.freq$count 
## [1] 1425857  766412   51335
  clas.sum3.freq$count 
## [1] 1540872  651397   51335
  clas.sum4.freq$count
## [1] 1580684  611585   51335
# Concordance frequencies decline as maps are layered upon one another.

Question #4

#* Save your data as R objects: # * All ensemble prediction maps as .img format #* Save these R objects in a .RData file as well

setwd(path.ex)
## Error in setwd(path.ex): object 'path.ex' not found
save(pers.BRT,pers.BRT2, file = 'ex12.RData')
save(pers.class.BRT, file = 'persClassBRT.img')
## Error in save(pers.class.BRT, file = "persClassBRT.img"): object 'pers.class.BRT' not found
save(pers.prob.BRT, file = 'persProbBRT.img')
## Error in save(pers.prob.BRT, file = "persProbBRT.img"): object 'pers.prob.BRT' not found

The End